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Abstract 

We  investigate  the  effects  of  smoothness  of  basis  functions  on  solution  accuracy  within 
the  isogeometric  analysis  framework.  We  consider  two  simple  one-dimensional  structural 
eigenvalue  problems  and  two  static  shell  boundary  value  problems  modeled  with  trivari¬ 
ate  NURBS  solids.  We  also  develop  a  local  refinement  strategy  that  we  utilize  in  one  of 
the  shell  analyses.  We  find  that  increased  smoothness,  that  is,  the  “k- method,”  leads  to  a 
significant  increase  in  accuracy  for  the  problems  of  structural  vibrations  over  the  classi¬ 
cal  C°-continuous  “p-method,”  whereas  a  judicious  insertion  of  Cf0-continuous  surfaces 
about  singularities  in  a  mesh  otherwise  generated  by  the  fc-method,  usually  outperforms  a 
mesh  in  which  all  basis  functions  attain  their  maximum  level  of  smoothness.  We  conclude 
that  the  potential  for  the  fc-method  is  high,  but  smoothness  is  an  issue  that  is  not  well  un¬ 
derstood  due  to  the  historical  dominance  of  (70-continuous  finite  elements  and  therefore 
further  studies  are  warranted. 

Key  words:  isogeometric  analysis,  finite  element  analysis,  ^-method,  p-method, 
refinement,  continuity,  smoothness,  structural  eigenvalue  problems,  shells,  singularities 
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1  Introduction 


The  concept  of  Isogeometric  Analysis,  introduced  by  Hughes,  Cottrell  and  Bazilevs 
[7]  and  further  developed  by  Cottrell  et  al.  [3]  and  Bazilevs  et  al.  [1],  was  initially 
motivated  by  the  gap  existing  between  Computer  Aided  Design  (CAD)  and  Finite 
Element  Analysis  (FEA).  The  first  manifestation  of  the  gap  is  in  the  initial  mesh 
generation  process.  The  design  is  encapsulated  in  some  type  of  CAD  model.  This 
model  often  includes  ambiguities,  such  as  gaps  and  overlaps,  and  levels  of  detail, 
such  as  individual  bolts,  welds,  etc.,  that  make  it  inappropriate  for  analysis.  The 
ambiguities  must  be  removed  and  defeaturing  must  be  performed  to  arrive  at  an 
Analysis  Suitable  Geometry  (ASG)  that  exactly  represents  the  features  of  interest 
for  the  calculation  (see  Figure  1).  This  ASG  must  then  be  replaced  with  a  finite 
element  mesh,  usually  a  piecewise  polynomial  approximation  of  the  actual  geom¬ 
etry.  Creating  a  mesh  can  be  one  of  the  more  time  consuming  steps  in  the  analysis 
process. 

Though  initial  mesh  generation  can  be  a  significant  bottleneck,  additional  difficul¬ 
ties  are  encountered  during  refinement.  Frequently,  if  an  accurate  solution  is  to  be 
obtained  through  a  series  of  refinements,  the  quality  of  the  geometric  approxima- 
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Remove  ambiguities,  such 
as  gaps  and  overlaps,  remove 
and/or  add  features,  construct 
one-  and  two-dimensional 
manifold  definitions  and  define 
section  properties,  etc. 

Fig.  1 .  The  geometry  of  an  object  of  engineering  interest  is  initially  encapuslated  in  a  Com¬ 
puter  Aided  Design  (CAD)  package.  The  CAD  description  must  frequently  be  changed 
significantly  to  create  an  Analysis  Suitable  Geometry  (ASG). 

tion  must  be  simultaneously  improved  or  else  the  error  will  reach  a  plateau  from 
which  it  cannot  be  reduced.  If  such  geometric  refinement  is  to  take  place,  a  link 
must  be  established  between  the  ASG  and  the  refinement  routine.  This  link  usually 
does  not  exist  in  practice  (see  Figure  2a).  This  may  be  one  of  the  reasons  why  auto¬ 
matic  refinement  has  had  little  impact  in  industry  despite  the  great  promise  shown 
in  academic  research  studies. 

Isogeometric  analysis  is  a  methodology  for  addressing  these  problems.  The  idea  is 
to  have  one  and  only  one  representation  of  the  geometry  which  exactly  encapsulates 
the  ASG  and  is  more  faithful  to  the  initial  CAD  representation.  While  an  ASG  must 
still  be  constructed,  using  functions  and  technologies  of  the  sort  found  in  CAD 
packages  may  facilitate  the  development  of  links  between  the  design  and  analysis 
software.  More  importantly,  if  the  finite  element  mesh  were  to  exactly  encapsulate 
the  ASG,  refinement  to  any  level  could  take  place  completely  within  the  analysis 
framework.  The  need  for  reestablishing  the  link  with  an  external  description  of  the 
geometry  would  be  completely  obviated  as  the  mesh  would  be  the  exact  geometry, 
as  in  Figure  2b. 

Our  current  implementation  of  the  isogeometric  analysis  concept,  based  on  Non- 
Uniform  Rational  B-Splines  (NURBS),  accomplishes  this  last  task  in  almost  all 
situations.  The  geometric  flexibility  of  the  NURBS  basis  allows  for  the  exact  repre¬ 
sentation  of  a  much  larger  class  of  objects  than  standard  finite  element  technology. 
Most  notably,  all  conic  sections  can  be  represented  exactly.  At  this  point,  genera¬ 
tion  of  the  initial  mesh  can  still  be  a  time  consuming  process  but  once  it  has  been 
performed,  the  isogeometric  mesh  encapsulates  the  exact  geometry  and  may  be 
refined  to  any  level  without  ever  altering  this  geometry  in  any  way. 

Meshless  methods  do  seem  to  share  certain  features  with  the  isogeometric  ap¬ 
proach.  The  description  of  complicated  geometries  within  such  methods,  however, 
has  been  almost  entirely  ignored  in  the  literature.  Notable  exceptions  are  found  in 
the  papers  of  Subbarayan  and  colleagues  [11,  20]  and  the  recent  work  of  Simkins 
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The  mesh  is  an 
approximate  encapsu¬ 
lation  of  the  ASG 


refinement  ? 


(a) 


Isogeometric 

Analysis 


The  mesh  is  an 
exact  encapsu¬ 
lation  of  the  ASG 


refinement 


(b) 

Meshless 


refinement  ? 


(c) 

Fig.  2.  The  analysis  process,  a)  In  finite  element  analysis,  mesh  refinement  requires  inter¬ 
action  with  an  external  description  of  the  geometry  if  the  quality  of  the  geometric  approx¬ 
imation  is  to  be  improved.  The  lack  of  such  interaction  is  an  impediment  to  adaptive  mesh 
refinement  procedures,  b)  In  isogeometric  analysis,  the  mesh  is  the  exact  geometry  and  so 
refinement  can  take  place  completely  within  the  analysis  framework,  c)  The  literature  on 
meshless  methods  is  yet  to  present  a  comprehensive  view  how  complex  geometries  may  be 
represented  and  how  that  representation  interacts  with  the  process  of  refining  the  solution 
space. 

et  al.  [15].  While  meshless  methods  do  show  great  promise  in  certain  areas,  a  clear 
view  of  the  proper  way  in  which  to  define  a  geometry,  as  well  as  how  that  descrip¬ 
tion  affects  both  refinement  of  the  solution  space  and,  perhaps  more  importantly, 
numerical  integration  of  the  basis  functions,  is  yet  to  emerge;  see  Figure  2c. 

An  important  feature  of  the  NURBS-based  approach  to  isogeometric  analysis,  that 
was  not  one  of  its  initial  motivations,  is  the  ability  to  use  functions  of  higher  order 
and  higher  continuity.  Section  2  will  describe  the  construction  of  NURBS  basis 
functions  that  may  have  up  to  p  —  1  continuous  derivatives  across  element  bound¬ 
aries,  where  p  is  the  order  of  the  underlying  polynomial.  This  is  seen  in  Section 
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3.1  to  have  a  profound  effect  in  structural  vibration  problems.  The  NURBS  func¬ 
tions  of  higher  continuity  offer  a  much  more  compact  representation  of  the  vibra¬ 
tional  modes  of  structures  than  do  standard  finite  element  functions,  yielding  much 
greater  accuracy  per  degree  of  freedom,  even  at  the  same  polynomial  order.  In  Sec¬ 
tions  3.2  and  3.3  we  study  shells  modeled  as  trivariate  NURBS  solids.  We  explore 
local  refinement  and  control  of  continuity.  Results  indicate  that  in  regions  with  very 
large  gradients,  the  use  of  functions  with  reduced  continuity  leads  to  more  accurate 
results  on  a  per-degree-of-freedom  basis,  at  least  on  coarse  meshes.  These  observa¬ 
tions  lead  to  the  conclusion  that  local  control  of  the  continuity  of  the  basis  is  a  tool 
to  be  exploited  in  efficiently  representing  many  types  of  solutions.  In  Section  4  we 
draw  conclusions. 


2  Overview  of  the  Isogeometric  Analysis  Framework 


Our  current  implementation  of  the  isogeometric  analysis  concept  is  based  on  Non- 
Uniform  Rational  B-Splines  (NURBS).  This  section  will  present  an  in  depth  dis¬ 
cussion  of  the  NURBS  functions  and  their  usage  in  representing  various  geometries 
comprised  of  a  single  NURBS  patch.  The  myriad  of  refinement  options  encompass¬ 
ing  classical  h-  and  p-refinement,  as  well  as  the  new  ^-refinement  described  in  [7], 
will  be  discussed  in  detail.  Lastly,  we  will  describe  the  use  of  multiple  patches  and 
local  refinement  using  constraint  equations. 


2. 1  B-splines  and  NURBS 

2.1.1  Knot  Vectors 

NURBS  are  built  from  B-splines  and  so  a  discussion  of  B-splines  is  a  natural  start¬ 
ing  point  for  the  investigation  of  NURBS.  Unlike  in  standard  FEA,  the  B-spline 
parametric  space  is  local  to  “patches”  rather  than  elements.  Patches  play  the  role 
of  subdomains  within  which  element  types  and  material  models  are  assumed  to  be 
uniform.  However,  a  variety  of  refinement  options  may  exist  within  a  single  patch. 
More  about  this  later.  Many  simple  domains  can  be  represented  by  a  single  patch. 

Note  that  the  distinction  between  “elements”  and  “patches”  may  be  thought  of  in 
two  different  ways.  In  [8]  and  [9],  the  patches  themselves  are  referred  to  as  el¬ 
ements.  This  is  not  unreasonable  as  the  parametric  space  is  local  to  patches  and 
a  finite  element  code  must  include  a  loop  over  the  patches  during  assembly.  As 
mentioned  previously,  we  take  the  alternate  view  that  patches  are  subdomains  com¬ 
prised  of  many  elements,  namely  the  “knot  spans”.  This  latter  view  seems  more 
appropriate  as,  in  our  current  code,  numerical  quadrature  is  being  carried  out  at  the 
knot  span  level.  Furthermore,  in  the  case  of  B-splines,  the  functions  are  piecewise 
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polynomials  where  the  different  “pieces”  join  along  knot  lines.  In  this  way  the  func¬ 
tions  are  C°°  within  an  element.  Lastly,  surprisingly  complicated  domains  may  be 
described  by  a  single  patch  ( e.g .,  all  of  the  numerical  examples  in  [7]).  Describing 
such  domains  as  being  comprised  of  one  element  seems  unnatural. 


A  knot  vector  in  one  dimension  is  a  set  of  coordinates  in  the  parametric  space, 
written  H  =  {£1,  £2, ...,  £n+p+i},  where  ^  6  M  is  the  ith  knot ,  i  is  the  knot  index, 
i  —  1,2,  ...,n  +  p  +  1,  p  is  the  polynomial  order,  and  n  is  the  number  of  basis 
functions  which  comprise  the  B-spline.  The  knots  partition  the  parameter  space 
into  elements.  Element  boundaries  in  the  physical  space  are  simply  the  images  of 
knot  lines  under  the  B-Spline  mapping,  as  shown  in  Figure  3. 


Fig.  3.  The  parametric  space  is  local  to  “patches”  rather  than  elements.  The  knots  partition 
the  patch  into  elements. 

Knot  vectors  may  be  uniform  if  the  knots  are  equally  spaced  in  the  parametric  do¬ 
main,  or  if  they  are  unequally  spaced,  they  are  non-uniform.  Knot  values  may  be 
repeated,  that  is,  more  than  one  knot  may  take  on  the  same  value.  The  multiplici¬ 
ties  of  knot  values  have  important  implications  for  the  continuity  properties  of  the 
basis.  A  knot  vector  is  said  to  be  open  if  its  first  and  last  knots  appear  p+1  times. 
Open  knot  vectors  are  standard  in  the  CAD  literature.  In  one  dimension,  basis  func¬ 
tions  formed  from  open  knot  vectors  are  interpolator  at  the  ends  of  the  parametric 
space  interval,  [£l>  £n+p+i]>  and  in  multiple  dimensions  they  are  interpolator  at  the 
comers  of  patches,  but  they  are  not,  in  general,  interpolator  at  interior  knots.  This 
is  a  distinguishing  feature  between  knots  and  “nodes”  in  finite  element  analysis. 


2.1.2  Basis  functions 


B-spline  basis  functions  are  defined  recursively  starting  with  piecewise  constants 

(P  =  0)  : 

'  1  if  &  <  £  <  £+1, 

0  otherwise. 


K.o(0  = 


(l) 


For  p  =  1,2,3,...,  they  are  defined  by 


NitP(0  =  J  NitP- 1(0  +  ^+1  /  Ni+ltP_ !($)■  (2) 


6+P  -  6 


6 


i+p+l  6+1 
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The  results  of  applying  (1)  and  (2)  to  a  uniform  knot  vector  are  presented  in  Figure 
4.  For  B-spline  functions  with  p  —  0  and  1,  we  have  the  same  result  as  for  standard 
piecewise  constant  and  linear  finite  element  functions,  respectively.  Quadratic  B- 
spline  basis  functions,  however,  are  different  than  those  in  FEM.  Each  quadratic  B- 
spline  is  identical  but  shifted.  This  distinguishes  them  from  quadratic  finite  element 
functions,  which  are  different  for  internal  and  end  nodes.  The  homogeneous  nature 
of  the  basis  has  implications  for  the  quality  of  the  approximation  and  the  potential 
for  efficient  solution.  In  the  case  of  structural  vibrations,  where  the  heterogeneity 
of  finite  element  functions  leads  to  a  branching  of  the  spectrum  that  degrades  the 
accuracy  of  a  large  percentage  of  the  computed  frequencies,  the  homogeneity  of  B- 
spline  functions  leads  to  dramatic  improvements,  as  will  be  shown  later  in  Section 


Fig.  4.  Basis  functions  of  order  0, 1,  2  for  uniform  knot  vector  S  =  {0, 1, 2, 3, 4, ...}. 


Fig.  5.  Quadratic  basis  functions  for  open,  non-uniform  knot  vector 

3  =  {0,0, 0,1,  2, 3, 4, 4,  5,  5,  5}. 

For  an  open,  non-uniform  knot  vector  we  can  attain  much  richer  behavior.  An  ex¬ 
ample  is  presented  in  Figure  5.  Note  that  the  basis  functions  are  interpolatory  at  the 
ends  of  the  interval  and  also  at  £  =  4,  the  location  of  a  repeated  knot,  where  only 
C'°-continuity  is  attained.  Elsewhere,  the  functions  are  C 1 -continuous.  In  general, 
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basis  functions  of  order  p  have  p  —  rrii  continuous  derivatives  across  knot  where 
rrii  is  the  multiplicity  of  the  value  of  £;  in  the  knot  vector.  When  the  multiplicity  of 
a  knot  value  is  exactly  p,  the  basis  is  interpolatory  there.  When  the  multiplicity  is 
p  +  1,  the  basis  becomes  discontinuous  and  the  patch  is  effectively  split  into  two 
separate  patches. 

An  important  property  of  B-spline  basis  functions  is  that  they  constitute  a  partition 
of  unity,  that  is,  V£, 

n 

E^UO  =  1-  (3) 

2=1 

This  is  a  feature  they  share  with  finite  elements  and  meshless  methods.  Also  of  note 
is  that  the  support  of  each  Ntp  is  compact  and  contained  in  the  interval  [£;,  £i+P+i] . 
Lastly,  observe  that  each  basis  function  is  point-wise  non-negative  over  the  entire 
domain,  that  is,  A^p(£)  >  0,  V£.  This  means  that  all  of  the  entries  of  a  mass  matrix 
will  be  positive,  which  has  implications  for  developing  lumped  mass  schemes. 


2.1.3  B-spline  curves 

B-spline  curves  in  are  constructed  by  taking  a  linear  combination  of  B-spline 
basis  functions.  The  vector- valued  coefficients  of  the  basis  functions  are  referred  to 
as  control  points.  These  are  analogous  to  nodal  coordinates  in  finite  element  analy¬ 
sis  in  that  they  are  the  coefficients  of  the  basis  functions,  but  the  non-interpolatory 
nature  of  the  basis  does  not  lead  to  the  usual  interpretation  of  the  control  point  val¬ 
ues.  Piecewise  linear  interpolation  of  the  control  points  gives  the  so-called  control 
polygon.  Again  note  that,  in  general,  control  points  are  not  interpolated  by  B-spline 
curves.  Given  n  basis  functions,  Npp,  i  =  1,2,  ...,n,  and  corresponding  control 
points  Bi  £  M.d,i  =  1,2,  ...,n,  a  piecewise-polynomial  B-spline  curve  is  given  by 

C(0  =  ENiA0Bi.  (4) 

2=1 

The  example  shown  in  Figure  6  is  built  from  the  quadratic  basis  functions  consid¬ 
ered  in  Figure  5.  The  curve  is  interpolatory  at  the  first  and  last  control  points,  a 
general  feature  of  a  curve  built  from  an  open  knot  vector.  Note  that  it  is  also  in¬ 
terpolatory  at  the  sixth  control  point.  This  is  due  to  the  fact  that  the  multiplicity 
of  the  knot  at  £  =  4  is  equal  to  the  polynomial  order.  Note  also  that  the  curve  is 
tangent  to  the  control  polygon  at  the  first,  last,  and  sixth  control  points.  The  curve 
is  Cp~ 1  =  C1  -continuous  everywhere  except  at  the  location  of  the  repeated  knot, 
£  =  4,  where  it  is  Cp~2  =  C° -continuous. 

The  properties  of  B-spline  curves  follow  directly  from  the  properties  of  their  basis 
functions.  For  example,  B-spline  curves  have  continuous  derivatives  up  to  order 
p  —  1  in  the  absence  of  repeated  knots  or  control  points.  Repeating  a  knot  or  control 
point  k  times  decreases  the  number  of  continuous  derivatives  by  k. 
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(a)  Curve  and  control  points  (b)  Curve  and  mesh  denoted  by  knot  locations 

Fig.  6.  B-spline,  piecewise  quadratic  curve  in  M2.  a)  Control  point  locations  are  denoted  by 
•’s.  b)  The  knots,  which  define  a  mesh  by  partitioning  the  curve  into  elements,  are  denoted 
by  B’s.  Basis  functions  and  knot  vector  as  in  Figure  5. 

Affine  transformations  of  a  B-spline  curve  are  obtained  by  applying  the  transfor¬ 
mations  directly  to  the  control  points.  This  turns  out  to  be  the  essential  property  for 
satisfying  so-called  “patch  tests,”  as  discussed  in  [7].  This  property  is  referred  to  as 
affine  covariance. 


2.1.4  h-refinement:  Knot  insertion 

The  mechanism  for  implementing  //-refinement  is  knot  insertion. 4  Knots  may  be 
inserted  without  changing  a  curve  geometrically  or  parametrically.  Given  a  knot 
Vector  =  {£i,  •••>  Cn+p+l})  1  =  {Cl  =  Cl)  ^2)  •••>  Cn+m+p+1  =  Cn+p+l} 

be  an  extended  knot  vector  such  that  E  C  E.  The  new  n  +  m  basis  functions  are 
formed  as  before  by  applying  (1)  and  (2)  to  the  new  knot  vector  E.  The  new  n  +  rn 
control  points,  B  =  { Bh  B2, ...,  Bn+m}T ,  are  formed  from  the  original  control 
points,  B  =  {f?i,  B2, ...,  Bn}T,  by 

B  =  TPB  (5) 

where 

rrtO 

1ij 

and 

+1  =  1+^-^  |z±£±1-|±i  for  g  =  0,  1,  2,  ...,p  —  1  (7) 

Ci+g+i  —  sj+i 

4  Note  that  in  the  CAD  literature  “knot  insertion”  refers  to  inserting  a  single  knot  into  a 
knot  vector,  whereas  “knot  refinement”  refers  to  inserting  multiple  knots  simultaneously. 
Here,  we  make  no  distinction  and  use  “knot  insertion”  to  refer  to  both  cases.  For  an  algo¬ 
rithm  for  inserting  an  individual  knot,  see  [7]. 


i  Ci  c  if  [Cj)£j+ 1) 

0  otherwise 
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Knot  values  already  present  in  the  knot  vector  may  be  repeated  as  above  but,  as 
described  subsequently  in  Section  2.1.2,  the  continuity  of  the  basis  will  be  reduced. 
Continuity  of  the  curve  is  preserved  by  choosing  the  control  points  as  in  (5),  (6)  and 
(7). 

Figure  7  shows  the  case  of  a  global  refinement  of  the  curve  from  Figure  6.  Insertion 
of  new  knot  values  has  parallels  with  the  classical  ^-refinement  strategy  in  finite 
element  analysis  as  it  splits  existing  elements  into  smaller  ones.  Repeating  existing 
knot  values  to  decrease  the  continuity  of  basis  does  not  have  an  analogue  in  FEA. 
We  will  return  to  this  idea  later. 


(a) 


(b) 


0,0,0  0.5  1  1.5  2  2.5  3  3.5  4,4  4.5  5,5,5 

2  =  {0, 0, 0,  .5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4, 4.5, 5, 5, 5} 


Fig.  7.  Knot  insertion,  a)  New  control  points  are  computed  from  the  original  control  points 
using  (5).  b)  Each  element  has  been  split  by  inserting  a  new  knot  at  the  midpoint  of  each 
knot  span. 

2.1.5  p-refinement:  Order  elevation 

The  mechanism  for  implementing  p-refinement  is  order  elevation  5  .  As  its  name 
implies,  the  process  involves  raising  the  polynomial  order  of  the  basis  functions 
used  to  represent  the  geometry  (and  the  solution  space,  as  our  finite  element  imple¬ 
mentation  will  be  isoparametric).  Recalling  from  Section  2.1.1  that  the  basis  has 
p  —  mi  continuous  derivatives  across  element  boundaries,  it  is  clear  that,  when  p  is 
increased,  ra;  must  also  be  increased  if  we  are  to  preserve  the  discontinuities  in  the 
derivatives  of  our  original  curve.  During  order  elevation,  the  multiplicity  of  each 
existing  knot  value  is  increased  by  one,  but  no  new  knot  values  are  added.  As  with 
knot  insertion,  neither  the  geometry  nor  the  parameterization  are  changed. 

5  sometimes  also  called  “degree  elevation.” 
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Fig.  8.  Order  elevation,  a)  New  control  points  are  calculated  so  as  to  preserve  the  geome¬ 
try  and  parameterization,  b)  The  mesh  remains  unchanged  as  no  new  elements  have  been 
created.  Note  the  increased  multiplicity  of  internal  knots.  This  is  done  to  preserve  disconti¬ 
nuities  in  the  derivatives  of  the  curve. 
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The  process  for  order  elevation  begins  by  replicating  existing  knots  until  their  mul¬ 
tiplicity  is  equal  to  the  polynomial  order,  thus  effectively  subdividing  the  curve  into 
many  Bezier  curves  by  knot  insertion  (see  Rogers  [14]  or  Farin  [5]  for  a  discussion 
of  Bezier  curves;  we  may  think  of  them  as  one  element  B-spline  curves).  The  next 
step  is  to  elevate  the  order  of  the  polynomial  on  each  of  these  individual  segments. 
Lastly,  excess  knots  are  removed  to  combine  the  segments  into  one,  order-elevated, 
B-spline  curve.  Several  efficient  algorithms  exist  which  combine  the  steps  so  as  to 
minimize  the  computational  cost  of  the  process.  Details  are  omitted  for  the  sake  of 
brevity.  For  a  thorough  treatment,  see  Piegl  and  Tiller  [12]. 

Figure  8  shows  this  process  applied  to  the  curve  in  Figure  6.  The  multiplicities  of 
the  knots  have  been  increased  but  no  new  elements  created.  Note  that  the  locations 
of  control  points  for  these  order-elevated  curves  are  different  than  those  in  the  h- 
refinement  example  (cf.  Figure  7). 


2.1.6  k-refinement:  Higher  order  and  higher  continuity 

As  we  have  seen,  the  two  primitive  operations  for  B-splines  are  knot  insertion  and 
order  elevation.  Knot  insertion  is  similar  to  //-refinement,  but  for  it  to  be  a  per¬ 
fect  analogue,  each  new  knot  value  would  have  to  be  inserted  with  multiplicity 
mi  =  p  to  ensure  a  C°  basis  everywhere.  Similarly,  if  we  begin  with  a  mesh  where 
all  functions  are  already  C°  across  element  boundaries,  order  elevation  coincides 
exactly  with  our  traditional  notion  of  //-refinement.  Knot  insertion  and  order  eleva¬ 
tion,  however,  provide  us  with  more  possibilities  than  the  two  standard  notions  of 
refinement. 

As  mentioned  above,  we  can  insert  new  knot  values  with  multiplicities  of  1  to  de¬ 
fine  new  elements  across  whose  boundaries  functions  will  be  Cp~ 1 .  We  can  also 
repeat  existing  knot  values  to  lower  the  continuity  of  the  basis  across  existing  el¬ 
ement  boundaries.  This  makes  knot  insertion  a  more  flexible  process  than  simple 
//-refinement.  Similarly,  we  have  a  more  flexible  higher-order  refinement  as  well.  It 
stems  from  the  fact  that  the  processes  of  order  elevation  and  knot  insertion  do  not 
commute.  If  a  unique  knot  value,  £,  is  inserted  between  two  distinct  knot  values  in 
a  curve  of  order  p,  the  number  of  continuous  derivatives  of  the  basis  functions  at 
£  is  p  —  1.  As  described  above,  if  we  subsequently  elevate  to  some  higher  order, 
q,  the  multiplicity  of  every  distinct  knot  value  (including  the  knot  just  inserted)  is 
incremented  ( q  —  p)  times  so  that  discontinuities  in  the  pth  derivative  of  the  basis 
are  preserved.  That  is,  the  basis  still  has  only  p  —  1  continuous  derivatives  at  £,  al¬ 
though  the  order  is  now  q.  If  instead  we  elevated  the  order  of  the  original,  coarsest 
curve  to  q  and  only  then  inserted  the  unique  knot  value  £,  the  basis  would  have  q  —  1 
continuous  derivatives  at  £.  We  refer  to  this  latter  procedure  as  k-refinement.  It  has 
no  analogue  in  standard  finite  element  analysis  6  . 

6  This  notion  of  //-refinement  is  not  the  same  as  the  “^-convergence”  described  in  [8]  in 
which  the  position  of  the  knots  is  altered.  It  bears  more  in  common  with  the  "/.’-version  finite 
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The  concept  of  ^-refinement  is  important  because  isogeometric  analysis  is  funda¬ 
mentally  a  higher-order  approach.  While  linear  finite  elements  can  be  represented 
within  a  NURBS  context,  it  takes  quadratic-level  NURBS  to  represent  conic  sec¬ 
tions  -  one  of  the  key  features  of  the  method.  In  traditional  p-refinement  there 
is  a  very  inhomogeneous  structure  to  arrays  due  to  the  different  basis  functions 
associated  with  surface,  edge,  vertex  and  interior  nodes.  In  addition,  there  is  a  pro¬ 
liferation  in  the  number  of  nodes  because  (7° -continuity  is  maintained  in  the  re¬ 
finement  process.  In  k-refinement,  there  is  a  homogeneous  structure  within  patches 
and  growth  in  the  number  of  control  variables  is  limited.  Let  us  emphasize  that 
an  “element”  in  one  dimension  is  defined  as  the  span  between  two  distinct  knot 
values.  The  number  of  elements  in  a  curve  will  then  be  the  number  of  non-zero 
knot  spans  in  the  knot  vector  ( e.g .,  the  domain  associated  with  the  knot  vector 
£  =  {0, 0, 0, 1, 2, 3, 3, 4, 4, 4}  consists  of  four  elements). 

Consider  the  classical  p-refinement  process.  Assume  the  initial  domain  consists  of 
one  element  and  p  +  1  basis  functions  (assuming  an  open  knot  vector),  which  we 
then  refine  by  inserting  new  knot  values  until  we  have  n  —  p  elements  and  n  basis 
functions,  all  Cp~l.  We  then  perform  order  elevation,  maintaining  continuity  at  the 
p—  1  level.  This  requires  replicating  each  distinct  knot  value,  adding  a  basis  function 
in  each  element  and  so  increasing  the  total  number  of  basis  functions  by  n  —  p  to 
2 n  —  p.  After  a  total  of  r  order  elevations  of  this  type,  we  have  (r  +  \)n  —  rp 
basis  functions,  where  p  is  still  the  order  of  our  original  basis  functions.  This  is 
seen  to  be  a  large  number  of  functions  when  one  considers  that  in  most  cases  of 
practical  interest  the  number  of  elements  will  be  quite  a  bit  larger  than  the  order  of 
the  basis.  By  comparison,  consider  beginning  with  the  same  one  element  domain 
and  proceed  by  k-refinement.  That  is,  order  elevate  r  times  adding  only  one  basis 
function  at  each  refinement,  then  insert  knots  until  we  have  n—p  elements  as  before. 
The  final  number  of  basis  functions  is  n  +  r,  each  having  r  +  p  —  1  continuity.  This 
amounts  to  an  enormous  savings  as  n + r  is  considerably  smaller  than  (r  + 1  )n — rp. 
Bear  in  mind  that  in  d  dimensions  these  numbers  are  raised  to  the  d  power.  Recall 
that  the  mesh,  defined  by  the  knot  locations ,  is  fixed  and  is  the  same  for  p-  and 
k-refinements.  See  Figures  9  and  10. 

It  is  important  to  note  that  “pure”  k-refinement,  where  all  functions  maintain  Cp~l 
continuity  across  element  boundaries,  is  only  possible  if  the  coarsest  mesh  is  com¬ 
prised  of  one  element.  If  the  initial  mesh  places  constraints  on  the  continuity  across 
certain  element  boundaries,  these  constraints  will  exist  on  all  meshes.  In  general, 
though  some  such  constraints  will  exist,  the  number  of  elements  desired  for  analysis 


element  method”  of  [16,  17]  in  that  k  refers  to  continuity,  but  the  motivations  are  different. 
The  increased  continuity  in  [16]  is  required  so  that  a  least-squares  finite  element  approach 
is  possible.  Such  an  approach  requires  that  the  solution  space  have  the  same  number  of 
continuous  derivatives  as  found  in  the  highest  order  derivative  of  the  differential  operator. 
Our  motivations  for  using  basis  functions  of  higher  continuity  are  efficiency  and  robustness 
of  the  solution  space  in  a  classical  Galerkin  finite  element  formulation  of  the  problem. 
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will  be  much  higher  than  the  number  needed  for  modeling  the  geometry.  Refine¬ 
ments  may  be  performed  such  that  the  functions  have  p  —  1  continuous  derivatives 
across  these  new  element  boundaries  and  the  benefits  of  ^-refinement  will  still  be 
significant. 


£  =  {0,0,l,l},p=l 


Fig.  9.  When  refining  a  coarse,  low-order  mesh  to  create  a  fine,  higher-order  mesh,  one 
may  choose  between  a  p-  or  fc-refinement  strategy.  Here  we  see  the  initial  step  for  each 
case,  (a)  Base  case  of  one  linear  element,  (b)  Classic  ^-refinement  approach:  knot  insertion 
is  performed  first  to  create  many  low-order  elements.  Subsequent  order  elevation  will  pre¬ 
serve  the  (7°  continuity  across  element  boundaries,  c)  New  ^-refinement  approach:  order 
elevation  is  performed  on  the  coarsest  discretization.  Subsequent  knot  insertion  will  result 
a  basis  which  is  Cp~l  across  the  newly  created  element  boundaries.  See  the  results  of  p- 
and  /e-refinement  for  several  different  polynomial  orders  in  Figure  10. 


2.1.7  The  hpk-refinement  space 

As  we  have  shown,  knot  insertion  and  order  elevation  are  the  primitive  operations 
by  which  classical  h-  and  ^-refinements,  as  well  as  the  new  k -refinement,  can  be 
implemented.  Recognizing  their  flexibility  as  compared  with  classical  refinement 
procedures  makes  feasible  the  notion  of  an  /ip/c-refinement  space.  Recalling  that  B- 
spline  curves  may  have  no  more  than  p—  1  continuous  derivatives  across  an  element 
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S  =  {0,0,0 ,1111,1,1,1}, p  =  2 


E  =  {0,0,  0,0,  §,  §,  §, 

|,  |,  |,  1, 1, 1, 1},  p  =  3 


=  10  0000  -  -  -  - 

u,  VJ,  u,  VJ,  3  ?  35  35 


22221111  ij  »=4 


H  =  {0,0,0,0,0,|, 

|,  1,1,1, 1,1},  p  =  4 


=  iO  00000----- 

w,  u,  u,  w,  u,  O)  O’)  O’)  35  O’) 


,  1,1, 1,1, 1,1},  p  =  5 


3  ’  3  ’  3  ’  3  ’  3 


(a) 


H  =  {0,0,0,0,0,0,|, 

§,l,l,l,l,l,l},p  =  5 

(b) 


Fig.  10.  Three  element,  higher-order  meshes  for  p-  and  /c-refinement.  a)  The  ^-refinement 
approach  results  in  many  functions  that  are  C°  across  element  boundaries,  b)  In  compar¬ 
ison,  /c-refinement  results  in  a  much  smaller  number  of  functions,  each  of  which  is  Cp~l 
across  element  boundaries. 
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boundary,  the  set  of  possible  refinements  may  be  characterized  as  in  Figure  11.  Pure 
^-refinement  keeps  h  fixed  but  increases  the  continuity  along  with  the  polynomial 
order,  as  in  Figure  12.  Pure  p-refinement  increases  the  polynomial  order  while  the 
basis  remains  C°,  as  in  Figure  13.  Increasing  the  multiplicity  of  existing  knot  val¬ 
ues  decreases  the  continuity  without  introducing  new  elements,  as  in  Figure  14. 
Inserting  new  knot  values  with  a  multiplicity  of  p  results  in  classical  /i-refinement, 
whereby  new  elements  are  introduced  that  have  C°  boundaries,  shown  in  Figure  15. 
Inserting  new  knot  values  with  a  multiplicity  of  1  decreases  h  without  decreasing 
the  minimum  continuity  already  found  in  the  mesh,  as  in  Figure  16.  Considering 
all  of  the  aforementioned  techniques  results  in  a  multitude  of  refinement  options 
beyond  simple  h-,  p-  and  ^-refinement,  see  Figure  17. 


Fig.  1 1 .  The  hpk- space.  The  set  of  all  allowable  refinements  is  contained  in  the  region 
shown  in  green.  Note  that  this  region  extends  in  the  direction  of  the  arrows. 


Fig.  12.  The  hpk- space.  In  pure  /.'-refinement,  the  locations  of  the  element  boundaries  (and 
thus  element  size,  h )  are  fixed.  As  the  polynomial  order,  p,  is  increased,  the  continuity  of 
the  functions  across  element  boundaries,  k,  is  increased  such  that  k  =  p  —  1  at  all  levels  of 
refinement. 
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Fig.  13.  The  hpk- space.  In  pure  p-refinement,  the  locations  of  the  element  boundaries  (and 
thus  element  size,  h)  are  fixed.  As  the  polynomial  order,  p,  is  increased,  the  continuity  of 
the  functions  across  element  boundaries,  k,  is  fixed  at  k  =  0  for  all  levels  of  refinement. 
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k  Knot  values  are  repeated 
4  -  to  reduce  A;  1 

3  -  1 


i 


Fig.  14.  The  hpk- space.  Repetition  of  existing  knot  values  decreases  the  continuity  across 
the  corresponding  element  boundary  without  creating  new  elements  or  changing  the  poly¬ 
nomial  order.  The  basis  has  p  —  m*  continuous  derivatives  across  knot  £*,  where  in,  is  the 
multiplicity  of  that  knot  value. 


Pure  /z-refinement 


Fig.  15.  The  hpk- space.  If  we  insert  new  knot  values  with  multiplicity  of  p,  new  elements 
are  created  and  the  basis  remains  C°  across  all  element  boundaries.  In  this  way  classical 
/j-refinement  is  exactly  replicated. 


2.1.8  Rational  B- splines 

As  described  in  the  beginning  of  this  section,  NURBS  are  formed  from  B-splines. 
Specifically,  NURBS  entities  in  Wl  can  be  obtained  by  projective  transformations 
of  B-spline  entities  in  M'i+I ,  in  particular,  conic  sections,  such  as  circles  and  el¬ 
lipses,  can  be  exactly  constructed  by  projective  transformations  of  piecewise  ra¬ 
tional  quadratic  curves.  The  projective  transformation  of  a  B-spline  curve  yields 
a  rational  polynomial  of  the  form  Cr{£)  =  /(£)/#(£) >  where  /  and  g  are  piece- 
wise  polynomials.  The  construction  of  a  rational  B-spline  curve  in  Md  proceeds 
as  follows.  Let  { B-v }  be  a  set  of  control  points  for  a  B-spline  curve  in  Md+ 1  with 
knot  vector  S.  These  are  referred  to  as  the  “projective  control  points”  for  the  de¬ 
sired  NURBS  curve  in  Wl.  The  control  points  in  Wl  are  derived  from  the  projective 


18 


Insertion  of  new  knot  values 


Fig.  16.  The  hpk- space.  Insertion  of  new  knot  values  with  a  multiplicity  of  1  results  in  a 
splitting  of  elements,  and  thus  a  decrease  in  h  (shown  in  the  figure  as  an  increase  in  /i_1). 
The  basis  has  p—  1  continuous  derivatives  across  these  new  element  boundaries,  and  so  the 
(possibly  lower)  minimum  continuity  already  existing  in  the  mesh  is  unchanged,  as  is  the 
polynomial  order. 


Fig.  17.  The  hpk- space.  Combining  knot  insertion  and  order  elevation  in  various  permuta¬ 
tions  allows  us  to  traverse  the  entire  allowable  refinement  space. 
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control  points  by  the  following  relations: 


(8) 

(9) 


w,  =  (BfW 1 


where  is  the  component  of  the  vector  etc.  and  w*  is  referred  to  as  the 


nth 


weight.  The  rational  basis  functions  and  NURBS  curve  are  given  by 
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Rational  surfaces  and  solids  are  defined  analogously  in  terms  of  the  rational  basis 
functions 
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The  powerful  thing  about  the  construction  of  the  NURBS  basis  functions  is  that, 
as  NURBS  in  are  B-splines  in  Md+1,  all  of  the  refinement  techniques  we  have 
discussed  are  applied  to  NURBS  by  operating  directly  on  those  higher  dimensional 
B-splines.  The  NURBS  basis  functions  also  form  a  partition  of  unity.  The  conti¬ 
nuity  and  support  of  NURBS  are  the  same  as  for  B-splines.  Affine  transformations 
in  physical  space  are  still  obtained  by  applying  the  transformation  to  the  control 
points,  that  is,  NURBS  possess  the  property  of  affine  covariance. 


2.2  Multiple  patches  and  local  refinement 


In  almost  all  practical  circumstances,  it  will  be  required  to  describe  a  domain  with 
multiple  NURBS  patches.  For  example,  if  different  material  or  physical  models  are 
to  be  used  in  different  parts  of  the  domain,  it  might  simplify  things  to  describe  these 
subdomains  by  different  patches.  Also,  if  different  subdomains  are  to  be  assembled 
in  parallel  on  a  multiple  processor  machine,  it  is  convenient  from  the  point  of  view 
of  data  structures  to  not  have  a  single  patch  split  between  different  processors. 
Most  common  is  the  case  where  the  domain  simply  differs  topologically  from  a 
cube.  The  tensor  product  structure  of  the  parameter  space  of  a  patch  makes  it  poorly 
suited  for  representing  complex,  multiply  connected  domains.  Such  geometries  can 
frequently  be  handled  quite  simply  by  using  multiple  patches  (see,  e.g..  Figure  18). 


20 


Fig.  18.  The  bracket  on  the  top  is  exactly  and  concisely  represented  by  five  simple  NURBS 
patches  (patch  boundaries  are  shown  in  red,  element  boundaries  in  blue).  The  patches  match 
geometrically  and  parametrically  on  the  internal  faces  where  they  meet. 
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Even  in  cases  where  a  cube  can  be  mapped  into  the  desired  object,  doing  so  might 
introduce  such  extreme  mesh  distortion  and  widely  varying  Jacobians  within  el¬ 
ements  that  analysis  will  be  adversely  affected.  Figure  19b  (from  [7])  shows  the 
amount  of  mesh  distortion  needed  to  represent  the  “stiffened  shell”  of  Figure  19a 
with  a  single  NURBS  patch.  A  mesh  using  multiple  patches,  shown  in  Figure  19c, 
exhibits  far  less  distortion  and  yields  a  much  more  “natural”  mesh. 


(a) 


(b)  (c) 

Fig.  19.  Multiple  patches  usually  produce  better  quality  meshes,  (a)  The  stiffened  shell  of 
[7]  can  be  modeled  using  a  single  NURBS  patch,  (b)  Such  a  mapping  produces  severe 
mesh  distortion  that  is  unavoidable  when  using  a  single  patch,  (c)  Allowing  the  shell  and 
the  stiffener  to  be  modeled  by  different  patches  creates  a  much  more  natural  mesh.  Patch 
boundaries  shown  in  red. 

Another  reason  for  using  multiple  patches  is  that  it  makes  local  refinement  possible. 
The  situation  is  represented  in  Figure  20.  Even  with  multiple  patches,  if  we  want 
the  control  points  of  the  two  patches  on  their  interface  to  be  in  one-to-one  corre¬ 
spondence,  we  need  to  have  matching  knot  vectors.  This  means  that  refinements  of 
one  patch  must  necessarily  propagate  from  that  patch  to  the  next.  If  we  are  to  allow 
knots  to  be  inserted  on  one  side  and  not  the  other  (/.<?.,  local  refinement),  we  may 
proceed  as  follows. 
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(a) 


(b) 


Fig.  20.  (a)  Global  refinement  employing  the  continuous  Galerkin  method,  (b)  Local  re¬ 
finement  employing  the  discontinuous  Galerkin  method  or  constraint  equations  at  the  patch 
level.  With  constraint  equations,  at  least  C° -continuity  can  be  attained  across  patches,  and 
higher-order  continuity  can  be  achieved  in  certain  cases  if  desired. 


Fig.  21.  The  two  patches  share  a  common  interface.  On  the  coarsest  mesh,  their  control 
points  on  that  interface  are  in  one-to-one  correspondence,  trivially  enforcing  C°  continuity. 


Fig.  22.  As  Patch  2  is  refined  by  knot  insertion  and  the  one-to-one  correspondence  of  the 
interface  control  points  is  lost.  Constraint  equations  may  be  employed  to  ensure  that  conti¬ 
nuity  is  maintained. 

Consider  the  two  B-spline  7  patches  that  meet  on  an  interface,  as  shown  in  Figure 
21.  On  the  coarsest  mesh,  we  assume  that  the  control  points  and  knot  vectors  in  the 
plane  of  the  face  are  identical  on  both  patches,  thus  ensuring  that  the  patches  match 
geometrically  and  parametrically  on  that  shared  face.  Using  superscripts  1  and  2  to 
identify  the  patch  numbers,  a  subscript  /  to  denote  control  points  on  the  face  where 
the  patches  meet,  and  a  subscript  n  to  denote  control  points  not  on  that  face,  we 


7  We  will  discuss  the  B-spline  case  here,  but  it  is  cmcial  to  note  that  if  we  were  to  use 
NURBS  rather  than  B-splines,  all  of  the  relationships  in  this  section  must  hold  for  the 
projective  control  points  and  projective  control  variables. 
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may  write  the  control  points  for  Patches  1  and  2  as 


B1 


VB9 


and  B2 


respectively,  where 


(14) 


(15) 


If  we  now  refine  the  basis  of  Patch  2  by  knot  insertion,  then  we  have  the  following 
new  set  of  control  points  for  Patch  2: 


B2  =  T  B2 


^Tn  0\ 

V°  T/J 


(16) 


where  T  is  the  multi-dimensional  generalization  of  the  extension  operator  defined 
in  (7).  As  before,  it  is  sparse  and  its  values  are  entirely  defined  by  the  knot  vectors 
and  the  polynomial  orders.  The  block  diagonal  structure  follows  from  the  fact  that 
we  are  using  open  knot  vectors.  When  open  knot  vectors  are  used,  each  face  of  a 
NURBS  solid  is  influenced  only  by  the  control  points  on  that  face.  Put  simply,  each 
face  of  the  NURBS  solid  is  a  NURBS  surface. 


Combining  (15)  and  (16),  we  see  that  C'°-continuity  of  the  geometry  is  maintained 
by  the  relationship 

B2f  =  TfB}.  (17) 

Building  on  the  approach  of  Kagan,  Fischer  and  Bar-Yoseph  [9] 8 ,  it  follows  that  for 
our  solution  space  to  enforce  the  same  continuity  constraints,  we  need  our  control 
variables  to  obey  precisely  the  same  relationship.  Let 


u 


l 


W/ 


and  u2 


(18) 


be  the  control  variables  on  Patch  1  and  the  refined  Patch  2,  respectively.  Then  C°- 
continuity  of  the  solution  across  the  interface  between  the  patches  may  be  main¬ 
tained  by  enforcing  the  constraint 

n)  =  Tfu).  (19) 


From  an  implementational  point  of  view,  the  two  patches  may  be  assembled  locally 
to  create  the  two  local  problems 


KW  =  b1  (20) 


8  In  [9],  a  similar  approach  was  taken  for  B -Splines  surfaces.  Here  we  extend  that  to 
NURBS  solids. 
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and 


K2u2  =  b2 


(21) 


for  the  control  points  on  either  patch.  Consistent  with  the  partitioning  of  the  control 
variables  in  (18),  we  partition  the  stiffness  matrices  as 


Kn  KiA 

K/n  K }f) 


Kn  Kf\ 

Kn  Kf) 


(22) 


Before  solving,  we  must  assemble  problems  (20)  and  (21)  into  one  global  problem 
accounting  for  the  behavior  of  both  patches,  as  well  as  their  interaction.  We  should 
have  three  coupled  blocks  of  equations:  one  corresponding  to  weighting  functions 
with  support  in  Patch  1  that  vanish  on  the  face  shared  by  the  two  patches,  one 
corresponding  to  weighting  functions  with  support  on  either  or  both  patches  that 
do  not  vanish  on  the  shared  face,  and  one  corresponding  to  weighting  functions 
with  support  on  Patch  2  that  vanish  on  the  shared  face.  We  begin  by  expanding  (20) 
using  the  partitioning  of  (22)  to  get 

Kn<  +  KKf  =  hn  (23) 

and 

K}„ui  +  K  ),n)  =  b‘.  (24) 

Inserting  (19)  into  (21)  and  expanding  yields 

K2nu2+K2/T/u}  =  b2  (25) 

and 

K2nu2+Kj/T/u}  =  b2.  (26) 

Note  that  (23)  is  the  block  of  equations  corresponding  to  weighting  functions  in 
Patch  1  that  vanish  on  the  shared  face.  Similarly,  (25)  is  the  block  of  equations  cor¬ 
responding  to  weighting  functions  in  Patch  2  that  vanish  on  the  shared  face.  Now 
(24)  and  (26)  both  correspond  to  weighting  functions  with  support  on  the  shared 
face  and  as  such  we  would  like  to  add  them  together  to  get  a  final  expression  for  that 
block.  Unfortunately,  they  contain  different  numbers  of  equations.  This  is  because 
we  assembled  the  two  patches  independently.  We  correctly  generated  the  equations 
in  (24)  by  testing  against  functions  in  the  “master”  weighting  space  associated  with 
Patch  1 ,  but  we  generated  the  equations  in  (26)  by  testing  against  all  of  the  functions 
in  the  larger  “slave”  weighting  space  on  Patch  2  without  regard  for  the  constraint. 
Just  as  the  basis  functions  of  the  slave  solution  space  on  Patch  2  corresponding  to 
the  shared  face  are  restricted  to  act  only  in  the  linear  combinations  defined  by  T/ 
that  result  in  functions  existing  in  the  master  solution  space,  so  too  must  the  func¬ 
tions  in  the  slave  weighting  space  act  only  in  such  linear  combinations  as  replicate 
functions  in  the  master  weighting  space.  This  constraint  may  be  enforced  by  now 
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premultiplying  (26)  by  Tj,  thus  constraining  the  weighting  functions  and  reducing 
the  number  of  equations  to  match  that  of  (24): 


TtK2  n2 

Lft^fnUn 


TtK2 


■/u/ 


■’tK2 


(27) 


We  may  now  express  the  global  system  comprised  of  (23),  (25),  and  ((24)+(27)j 


as 


where 


and 


K  = 


Ku  =  b, 


Kf 


K/n  (K}y  +  T }K}fT j)  TjK2n 


L fn 

0 


u 


Kn/T/ 


u1 


ui 


Vu"7 


K2 


b  = 


b1 


b} 


fTjb2 

b2 


(28) 


(29) 


(30) 


(31) 


We  may  recover  m  via  (19)  after  solving  (28). 


This  approach  ensures  C°  continuity  in  the  solution  across  the  patch  boundary  when 
one  patch  is  a  knot  refined  version  of  the  other  patch  on  their  common  interface. 
Higher  continuity  has  also  been  implemented  by  applying  similar  constraint  equa¬ 
tions  in  the  normal  direction.  As  long  as  the  geometries  are  compatible,  the  patch 
boundary  may  be  seen  as  the  result  of  inserting  a  knot  into  some  “metapatch”  P+1 
times.  It  should  be  noted  that  these  are  strong,  exact  constraints,  not  approxima¬ 
tions.  An  approach  that  would  allow  for  weak  enforcement  of  continuity,  as  well 
as  allowing  for  local  order  elevation  is  to  use  discontinuous  Galerkin  techniques 
at  the  patch  level.  That  is,  weakly  enforce  continuity  of  appropriate  fluxes  across 
patch  boundaries  while  strongly  enforcing  them  across  element  boundaries  within 
the  patch. 


Remark 


It  is  important  to  note  that  these  operations  could  also  be  applied  over  the  entire 
domain  rather  than  just  for  the  interface  between  patches.  These  could  be  used 
in  a  multigrid  scheme  where  the  grid  transfer  operator  would  be  T.  This  could 
potentially  be  very  efficient  as  T  is  uniquely  defined  by  the  knot  vectors  and  thus 
its  construction  is  very  inexpensive. 
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3  Numerical  examples 


In  Section  2. 1 .6  we  compared  the  number  of  degrees-of-freedom  in  k-  and  p-refined 
meshes.  We  found  that,  for  the  same  mesh  and  polynomial  order,  ^-refinement  in¬ 
volved  many  fewer  degrees-of-freedom  than  p-refinement.  This  suggests  to  us  that 
/c-refinement  may  be  a  more  efficient  procedure  than  p-refinement.  However,  this  is 
not  completely  clear  because  other  factors  are  at  play.  A  traditional  way  to  assess 
efficiency  is  by  comparing  accuracy  on  a  per  degree-of-freedom  basis,  although 
this  may  not  be  entirely  satisfactory  either.  Nevertheless,  to  get  some  sense  of  the 
relative  efficiency  of  /c -refined  and  p-refined  meshes,  we  will  adopt  this  approach  in 
the  following  numerical  examples.  We  will  often  refer  to  p-refined  meshes  simply 
as  “finite  elements”  and  /--refined  meshes  as  “NURBS.” 


3.1  Vibrations  of  beams  and  rods 


We  study  the  problem  of  the  structural  vibrations  of  an  elastic  fixed-fixed  rod  of  unit 
length,  whose  natural  frequencies  and  modes,  assuming  unit  material  parameters, 
are  governed  by: 

u,xx  +  u2u  =  Ofora:  e  ]0, 1[ 

u(0)  =  u(l)  =  0,  (  j 

and  for  which  the  exact  natural  frequencies  are: 

ujn  =  nn,  with  n  =  1,2,  3...  (33) 


As  a  first  numerical  experiment,  the  eigenproblem  is  solved  with  both  finite  el¬ 
ements  and  isogeometric  analysis  using  quadratic  basis  functions.  The  resulting 
natural  frequencies,  are  presented  in  Figure  23,  normalized  with  respect  to  the 
exact  solution  (33),  and  plotted  versus  the  mode  number,  n,  normalized  by  the  total 
number  of  degrees-of-freedom,  N.  To  produce  the  spectra  of  Figure  23,  we  used 
N  =  999  but  the  results  are  in  fact  independent  of  N. 

Figure  23  illustrates  the  superior  behavior  of  NURBS  basis  functions  compared 
with  finite  elements.  In  this  case,  the  finite  element  results  depict  an  acoustical 
branch  for  n/N  <0.5  and  an  optical  branch  for  n/N  >  0.5  (see  Brillouin  [2]).  As 
we  go  to  higher-order,  the  disparity  becomes  even  greater.  Higher-order  NURBS 
outperform  higher-order  finite  elements  by  an  ever  increasing  margin,  see  Figure 
24. 

Additionally,  transverse  vibrations  of  a  simply-supported,  unit  length  Bemoulli- 
Euler  beam  are  considered  (see  Hughes  [6],  Chapter  7).  For  this  case,  the  natural 
frequencies  and  modes,  assuming  unit  material  and  cross-sectional  parameters,  are 
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1.3 


Fig.  23.  Fixed-fixed-rod.  Normalized  discrete  spectra  for  quadratic  finite  elements  and 
NURBS. 


Fig.  24.  Fixed-fixed-rod.  Normalized  discrete  spectra  for  higher-order  finite  elements  and 
NURBS. 
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governed  by: 


^ ,xxxx  ic  u  0  for  x  £  ]0, 1[  ,3., 

u(  0)  =  u(  1)  =  utXX(  0)  =  ujXX(  1)  =  0, 

where 

uin  =  (mi)2,  with  n  =  1,  2,  3, ...  (35) 

The  numerical  experiments  and  results  for  the  Bemoulli-Euler  beam  problem  are 
analogous  to  the  ones  reported  for  the  rod.  Note  that  the  classical  beam  finite  el¬ 
ement  employed  to  solve  problem  (34)  is  a  two-node  Hermite  cubic  element  with 
two  degrees-of-freedom  per  node  (transverse  displacement  and  rotation),  whereas 
our  isogeometric  analysis  formulation  is  rotation-free  (see,  for  example,  Engel  et 
al.  [4]).  Figure  25  presents  the  discrete  spectra  obtained  using  different  order  finite 
element  and  NURBS  basis  functions.  Again,  ^-refinement  results  are  dramatically 
better  on  a  per  degree-of-freedom  basis. 


Fig.  25.  Simply-supported  beam.  Normalized  discrete  spectra  for  higher-order  finite  ele¬ 
ments  and  NURBS. 

Remark 

It  is  very  important  to  observe  the  trends  in  Figures  24  and  25.  For  finite  elements, 
the  optical  branches  of  the  frequency  spectra  diverge  as  p  is  increased.  That  is,  the 
errors  in  the  higher  frequencies  get  worse  as  p  is  increased.  It  is  well-known  that 
higher  frequencies  are  inaccurate  in  finite  element  analysis,  but  it  is  apparently  a 
new  observation  that  they  get  progressively  worse  as  p  is  increased.  On  the  other 
hand,  for  NURBS,  the  entire  spectrum  converges  as  p  is  increased.  These  oppo¬ 
site  trends  may  be  very  important  in  applications  such  as  wave  propagation  and 
turbulence,  in  which  the  entire  discrete  spectrum  may  participate  significantly  in 
the  solution.  We  conjecture  that  NURBS,  capable  of  attaining  almost  spectral  ac¬ 
curacy  on  patches,  as  evidenced  by  Figures  24  and  25,  may  be  superior  to  classical, 
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higher-order,  C'°-continuous  finite  elements  in  these  applications.  It  may  also  be 
noted  that,  based  on  similar  studies,  NURBS  exhibit  superior  accuracy  compared 
to  finite  elements  for  first-order  spatial  operators.  This  buttresses  the  belief  that 
NURBS  should  be  capable  of  attaining  better  accuracy  than  finite  elements  in  rep¬ 
resenting  wave  phenomena  and  turbulence. 


3.2  Hyperboloidal  shell 


The  hyperboloidal  shell  problem  was  introduced  to  us  by  Prof.  Bama  Szabo.  His 
group  had  analyzed  the  structure  using  a  /> refinement  strategy  on  meshes  created 
using  quasi-regional  mappings  based  on  optimal  collocation  at  Babuska  points  [18]. 
He  was  interested  in  an  independent  estimate  of  the  limit  value  of  the  potential  en¬ 
ergy  and  suggested  we  investigate  it  as  our  isogeometric  approach  is  capable  of 
exactly  representing  the  conic  sections  in  the  geometry.  The  problem  was  consid¬ 
ered  previously  by  Lee  and  Bathe  [10]  using  shell  elements,  but  unfortunately  they 
do  not  report  the  potential  energy  in  their  results. 

The  domain  is  the  thin-walled  solid  seen  in  Figure  26,  whose  mid-surface  is  defined 
by 

x2  +  z2  -y2  =  1,  ye  [-1,1].  (36) 

The  structure  has  a  thickness  of  t  —  0.001  in  the  direction  normal  to  this  mid¬ 
surface  (all  distances  are  in  meters).  The  loading  is  a  smoothly  varying  pressure 
normal  to  the  surface, 

p(6)  =p0  cos(2#),  (37) 

with  p0  =  1.0  MPa.  The  top  and  bottom  of  the  structure  are  fixed. 


X 


Fig.  26.  The  geometry  of  the  hyperboloidal  shell. 
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3.2.1  Mesh  generation  and  implementation 

Only  a  quarter  of  the  structure  is  modeled  due  to  symmetry.  The  mid-surface  is  a 
conic  section,  namely  a  hyperbola,  extruded  in  a  path  defined  by  another  conic  sec¬ 
tion,  a  circle.  As  rational  quadratic  NURBS  are  capable  of  representing  all  conic 
sections,  this  hyperboloidal  surface  of  revolution  can  be  represented  exactly.  How¬ 
ever,  the  inner  and  outer  surfaces  of  the  structure  are  defined  as  offsets  of  the 
mid-surface,  shifted  by  ±t/ 2  in  the  normal  direction,  and  are  not  conic  sections. 
Moreover,  they  are  not  in  the  NURBS  space,  so  our  mesh  will  inherently  be  an 
approximate  geometry  9 . 

The  decision  was  made  to  use  two  quadratic  elements  through  the  thickness  of 
the  structure  (see  Figure  27).  The  knot  value  defining  the  boundary  between  the 
elements  has  a  multiplicity  equal  to  its  polynomial  order,  2,  thereby  making  the 
geometrically  exact  mid-surface  a  discernible  entity  within  the  mesh  -  it  is  the 
boundary  between  the  inner  and  outer  layers  of  elements.  Knots  are  then  inserted 
into  the  appropriate  knot  vectors  to  define  the  elements  in  the  mid-surface  of  the 
coarsest  mesh.  This  mid-surface  mesh  is  identical  for  all  polynomial  orders. 


Fig.  27.  Two  quadratic  elements  are  used  through  the  thickness.  The  basis  is  C°  across 
the  interior  element  boundary,  thus  making  the  boundary  itself  a  NURBS  surface.  By  this 
construction,  the  geometrically  exact  mid-surface  is  a  discernible  entity  within  the  mesh. 

Once  the  coarse  mid-surface  mesh  is  fixed,  so  too  is  the  number  of  basis  functions 
in  the  axial  direction.  The  offset  curves  that  define  the  inner  and  outer  surfaces  of 
revolution  must  now  be  interpolated.  The  number  of  points  along  each  curve  that 
may  be  interpolated  is  equal  to  the  number  of  basis  functions  in  the  axial  direc¬ 
tion.  Due  to  the  use  of  open  knot  vectors  (see  [7]),  the  number  of  functions  for  a 
fixed  mesh  grows  with  p.  Specifically,  for  the  chosen  mesh  there  are  14  +  p  ba¬ 
sis  functions  in  the  axial  ( i.e .,  y)  direction.  As  a  result,  14  +  p  points,  equispaced 
in  the  parametric  domain,  are  calculated  along  the  hyperbola,  then  offset  by  ±t/ 2 
using  the  analytically  computed  normals  to  the  curve.  These  offset  points  are  then 
interpolated  using  C p~l  B-splines  to  create  the  approximate  geometry.  In  this  way, 
the  quality  of  the  overall  geometric  approximation  improves  as  the  polynomial  or- 


9  Note  that  it  is  the  offset  of  the  hyperbola  which  is  not  represented  exactly.  In  the  radial 
direction,  the  offset  of  a  circle  is  again  a  circle  and  therefor  exists  in  our  NURBS  space.  It 
is  the  radii  of  the  circles  denoting  the  inner  and  outer  surfaces  of  the  structure  at  a  given 
height  y  that  are  not  exact. 
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der  increases,  though  the  mid-surface  mesh  is  the  same  for  all  orders.  The  loading, 
however,  does  not  differ  as  it  is  applied  directly  to  the  mid- surface  itself. 

To  complete  Mesh  1  for  each  polynomial  order,  knots  are  inserted  near  the  fixed 
ends  creating  two  rows  of  small  elements  in  order  to  better  resolve  the  boundary 
layer.  The  multiplicity  of  these  knots  is  p,  and  so  the  basis  functions  are  C°  across 
these  element  boundaries,  shown  in  red  in  Figure  28.  As  discussed  above,  we  could 
have  introduced  the  same  number  of  new  degrees-of-freedom  into  the  region  by 
creating  many  small  elements,  each  having  p  —  1  continuous  derivatives  across 
their  boundaries.  Instead  we  have  chosen  fewer  elements  with  lower  continuity. 
The  motivation  for  introducing  these  C°  mesh-lines  is  that  previous  experience  has 
indicated  that  doing  so  helps  to  prevent  the  behavior  in  the  layer  from  polluting 
results  elsewhere  in  the  domain.  The  main  reason  for  this  is  that  introducing  a  C° 
mesh-line  results  in  a  localizing  of  the  support  of  the  basis  functions  and  thereby 
decreasing  the  coupling  of  functions  within  the  boundary  layer  with  functions  out¬ 
side  of  the  boundary  layer  (see  Figure  29).  The  result  is  crisper  layers  and  more 
compact  representations  of  the  global  solution,  particularly  on  coarse  meshes 10 . 

Meshes  2  and  3,  seen  in  Figures  30  and  31,  respectively,  are  the  result  of  subsequent 
uniform  fi-refinements.  The  basis  is  Cp~x  across  the  element  boundaries  introduced 
through  these  refinements.  The  geometry  and  parameterization  remain  unchanged, 
and  so  Mesh  1  fixes  the  geometry  for  each  polynomial  order.  In  this  way,  our  anal¬ 
ysis  is  an  fi-method,  repeated  for  several  different  polynomial  orders,  rather  than 
a  p-method  repeated  for  several  meshes  as  in  [18].  Recall,  however,  that  the  mid¬ 
surface  meshes  for  Mesh  1,  Mesh  2  and  Mesh  3  are  independent  of  the  polynomial 
order.  This  was  done  in  an  effort  to  make  the  results  from  one  polynomial  order  to 
the  next  as  comparable  as  possible. 


3.2.2  Results 

The  numbers  of  degrees-of-freedom  and  the  numerically  computed  volume  are  re¬ 
ported  in  Table  1 .  The  potential  energy  for  each  mesh,  as  well  as  the  limit  as  h  — »  0 
estimated  using  Richardson  extrapolation  (see  the  appendix),  are  reported  in  Table 
2.  Plots  of  the  deformed  geometry  as  seen  from  several  different  angles  are  shown 
in  Figures  32-34.  The  displacement  has  been  amplified  by  a  factor  of  10  to  make  it 
more  visible.  Due  to  the  sinusoidal  character  of  the  loading,  the  deformed  structure 
has  “compression  lobes”  and  “expansion  lobes.”  In  both  cases,  the  largest  gradients 
of  the  solution  are  contained  in  thin  layers  near  the  fixed  ends  of  the  structure.  Plots 

10  Note  that  we  would  expect  a  very  fine  mesh  (e.g.,  one  with  all  of  the  elements  the  size 
of  those  in  the  boundary  layer)  comprised  of  highly  continuous  functions  to  represent  the 
solution  more  efficiently  than  a  classical  p-method  on  a  per  degree-of-freedom  basis.  On  a 
coarse  mesh,  however,  the  efficiency  of  the  method  is  degraded  if  a  large  percentage  of  the 
functions  have  support  in  both  very  small  and  very  large  elements.  This  issue  is  investigated 
in  more  detail  in  Section  3.3. 
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Fig.  28.  Mesh  1.  Basis  functions  have  p  —  1  continuous  derivatives  across  blue  element 
boundaries.  They  are  only  C°  across  red  element  boundaries. 

of  the  radial  and  vertical  displacement  at  a  compression  lobe  and  at  an  expansion 
lobe  are  shown  in  Figures  35-42.  In  all  of  these  plots,  results  for  each  of  the  polyno¬ 
mial  orders  on  Mesh  3  are  shown.  While  the  quadratics  are  far  from  converged,  and 
cubics  seem  to  be  showing  signs  of  the  geometry  error,  the  quartics  and  quintics  lie 
practically  on  top  of  each  other.  It  is  likely  that  if  the  study  were  extended  to  higher 
polynomial  orders,  the  curves  would  be  virtually  indistinguishable  from  the  quintic 
solution. 


Mesh  1  DOF 

Mesh  2  DOF 

Mesh  3  DOF 

Vol.  of  shell 

2,2 

2160 

6300 

21060 

1.597535  *  10“2  m3 

3,2 

3045 

7755 

23655 

1.597527  *  10“2  m3 

4,2 

4080 

9360 

26400 

1.597530  *  10"2  m3 

5,2 

5265 

11115 

29295 

1.597530  *  10"2  m3 

e  1 

Mesh  data.  Here  p  denotes  the  polynomial  order  in  the  plane  of  the  surface,  while  q  is  the 
polynomial  order  through  the  thickness.  The  exact  volume  of  the  shell  is  1.597530  *  10-2 
m3. 

After  this  initial  study  was  completed,  a  second  study  was  performed  using  the 
maximum  continuity  possible.  The  shell  geometries  were  identical  to  those  pre¬ 
sented  above,  as  were  the  meshes  outside  of  the  boundary  layer  region.  The  width 
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E  =  {0,0, 0,0, 0.1, 0.1, 0.1, 1,1, 1,1} 

(a) 


E  =  (0, 0, 0, 0, 0.03, 0.06, 0.1, 1, 1, 1, 1} 

(b) 

Fig.  29.  Example  boundary  layer  meshes.  Both  meshes  have  the  same  number  of  basis  func¬ 
tions.  a)  Only  one  of  the  seven  basis  functions  has  support  both  inside  the  layer  (£  <  0.1) 
and  outside  of  the  layer  (£  >  0.1).  b)  Three  of  the  seven  basis  functions  have  support  both 
inside  and  outside  of  the  layer. 


Mesh  1 

Mesh  2 

Mesh  3 

Est.  limit 

2,2 

-4.668902 

-4.751145 

-4.796779 

-4.799821 

3,2 

-4.783082 

-4.794395 

-4.801991 

-4.802112 

4,2 

-4.787948 

-4.795994 

-4.799878 

-4.799893 

5,2 

-4.791334 

-4.798373 

-4.799941 

-4.799942 

*  10-2  MNm 

*  10"2  MNm 

*  10-2  MNm 

*  10"2  MNm 

Table  2 

Potential  energy.  Estimated  limit  calculated  using  Richardson  extrapolation. 
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Fig.  30.  Mesh  2.  The  second  mesh  is  generated  by  uniform  /z-refinement  of  Mesh  1.  The 
basis  is  Cv~x  across  the  new  element  boundaries. 


Fig.  31.  Mesh  3.  The  third  mesh  is  generated  by  uniform  /z-refinement  of  Mesh  2.  The  basis 
is  Cv~x  across  the  new  element  boundaries. 
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Fig.  32.  The  deformed  configuration.  Displacements  have  been  amplified  by  a  factor  of  10 
for  visibility. 


Fig.  33.  The  deformed  configuration.  Compression  lobes  are  visible  where  the  loading  is 
directed  inward. 


Fig.  34.  The  deformed  configuration.  Expansion  lobes  are  visible  where  the  loading  is  di¬ 
rected  outward. 
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Radial  Displacement  v.  Height  at  Compression  Lobe 


Fig.  35.  Compression  lobe.  Radial  displacement  versus  height. 


Radial  Displacement  v.  Height  at  Compression  Lobe 


y 

Fig.  36.  Compression  lobe.  Detail  of  radial  displacement  versus  height. 


Vertical  Displacement  v.  Height  at  Compression  Lobe 


Fig.  37.  Compression  lobe.  Vertical  displacement  versus  height. 
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x  10  3  Vertical  Displacement  v.  Height  at  Compression  Lobe 


Fig.  38.  Compression  lobe.  Detail  of  vertical  displacement  versus  height. 


Radial  Displacement  v.  Height  at  Compression  Lobe 


Fig.  39.  Expansion  lobe.  Radial  displacement  versus  height. 


Radial  Displacement  v.  Height  at  Compression  Lobe 


Fig.  40.  Expansion  lobe.  Detail  of  radial  displacement  versus  height. 
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Vertical  Displacement  v.  Height  at  Expansion  Lobe 


Fig.  41.  Expansion  lobe.  Vertical  displacement  versus  height. 


x  10 


-3  Vertical  Displacement  v.  Height  at  Expansion  Lobe 
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Fig.  42.  Expansion  lobe.  Detail  of  vertical  displacement  versus  height. 
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functions  may  result  in  a  more  accurate  solution. 


Mesh  1 

Mesh  2 

Mesh  3 

C°  Boundary  layer 

0.1804% 

0.0337% 

0.0011% 

Cp~l  Boundary  layer 

0.3379% 

0.0206% 

0.0007% 

Table  3 

Comparison  of  the  potential  energy  errors  in  the  two  approaches  to  boundary  layer  meshing, 
p  =  5.  Though  the  coarse  mesh  favors  the  C°  boundary  layer,  smooth  functions  prove  to 
be  more  accurate  once  the  meshes  are  sufficiently  fine.  The  reference  solution  used  in  the 
error  calculation  is  the  estimated  limit  for  the  quintic  case  from  the  previous  table. 

Remark 

In  preliminary  investigations  of  this  problem,  a  mesh  was  used  that  had  one  cubic 
element  through  the  thickness,  rather  than  the  two  quadratic  elements  as  described 
above.  The  interpolation,  and  thus  the  geometry  itself,  was  identical.  The  pressure 
loading  was  applied  to  the  inner  surface  of  the  structure.  Though  a  rigorous  study 
was  not  performed  under  these  conditions,  it  was  clear  that  the  potential  energy 
was  consistently  greater  (less  negative)  by  about  3%  than  the  limit  values  reported 
by  Szabo  [18].  The  meshes  used  in  the  present  study  were  generated  in  an  effort  to 
improve  upon  these  early  efforts,  as  indeed  they  did.  Pressurizing  the  mid-surface 
instead  of  the  inner  surface  not  only  allowed  the  use  of  accurate  normals  (which 
was  seen  as  the  most  likely  source  of  error  in  the  preliminary  efforts),  but  also 
made  the  study  more  consistent  as  the  loading  is  now  the  same  for  all  polynomial 
orders. 

When  comparing  with  the  results  of  Szabo  [18],  our  value  for  the  potential  energy 
is  slightly  more  negative  but  the  difference  is  only  ~  0.1%  of  the  total  energy.  It  is 
likely  that  the  discrepancy  is  due  to  geometrical  differences  in  our  models,  though 
the  implementations  of  the  loading  may  have  had  an  effect  as  well.  Despite  having 
the  exact  geometrical  mid-surface,  our  approach  to  interpolation  was  somewhat 
arbitrary,  though  the  volume  calculations  imply  that  our  geometries  for  p  =  4  and 
p  =  5  are  more  accurate  than  those  in  [18]. 

3.3  Hemispherical  shell  with  a  stiffener 


The  hemispherical  shell  with  a  stiffener  problem  (see  Figure  43)  was  modeled  with 
a  single  NURBS  patch  in  [7],  As  was  shown  in  Section  2.2,  use  of  a  single  patch 
leads  to  substantial  distortion  of  the  elements.  While  it  speaks  well  of  the  overall 
robustness  of  the  method  that  accurate  results  were  still  obtained,  efficiency  clearly 
suffered.  The  p-method  used  in  that  convergence  study  was  not  competitive  on 
a  per  degree-of-freedom  basis  with  the  original  results  of  Rank  et  al.  [13],  who 
used  a  trunk  space  p-refinement  strategy.  Such  an  approach  does  not  use  the  full 
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tensor  product  space  of  basis  functions,  but  the  much  smaller  trunk  space,  just  large 
enough  to  ensure  the  optimal  convergence  rate  at  a  given  polynomial  order  (see 
Szabo,  Duster  and  Rank  [19]  for  a  discussion  of  the  trunk  space  and  the  p-method 
in  general).  As  NURBS  necessarily  have  an  underlying  tensor  product  structure,  at 
least  on  patches,  an  analogous  isogeometric  analysis  approach  exploiting  the  trunk 
space  has  not  been  attempted  thus  far. 


Point  B 
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Fig.  43.  Hemispherical  shell  with  stiffener.  Problem  description  from  Rank  et  al.  [13]. 


Despite  the  tensor  product  structure  of  NURBS,  ^-refinement  presents  the  possi¬ 
bility  of  improved  efficiency.  In  fact,  ^-refinement,  in  conjunction  with  the  use  of 
multiple  patches  to  create  better  quality  meshes,  and  the  use  of  local  refinement  to 
avoid  placing  functions  in  regions  where  they  are  not  needed,  enables  the  NURBS 
based  approach  to  show  an  accuracy  per  degree-of-freedom  comparable  to  the  re¬ 
sults  presented  by  Rank  et  al.  in  [13];  see  Figures  46-49. 


In  studying  the  stiffened  shell  problem  with  a  ^-refinement  approach,  certain  pre¬ 
viously  unobserved  features  begin  to  emerge.  As  above,  if  a  given  mesh  of  higher- 
order  and  high  continuity  does  not  achieve  the  level  of  accuracy  desired,  one  can 
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Fig.  44.  Hemispherical  shell  with  stiffener.  The  coarse  mesh  may  be  refined  in  multiple 
ways. 
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(a)  Local  /^-refinement 


(b)  Local  k* -refinement 


Fig.  45.  Hemispherical  shell  with  stiffener,  (a)  A  ^-refinement  approach  with  Cv~x  conti¬ 
nuity  across  element  boundaries.  Many  small  elements  are  used  to  get  a  well-resolved  so¬ 
lution.  (b)  Functions  are  C°  across  the  element  boundaries  in  red,  Cp~l  elsewhere.  Fewer 
elements  are  needed  than  in  (a).  In  both  cases,  the  basis  is  C°  across  patch  boundaries, 
shown  in  black,  and  local  refinement  is  implemented  at  the  patch  level. 
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add  more  degrees-of-freedom  by  inserting  a  knot  in  one  of  the  parametric  direc¬ 
tions.  The  number  of  new  degrees-of-freedom  is  exactly  the  same  regardless  of 
whether  a  new  knot  value  is  inserted  (creating  new  elements  by  splitting  existing 
ones),  or  whether  an  existing  knot  value  is  repeated  (creating  no  new  elements, 
but  decreasing  the  continuity  of  the  basis  across  the  corresponding  element  bound¬ 
aries).  While  a  rigorous  analysis  of  the  two  approaches  has  not  yet  been  performed, 
in  the  present  results  it  seems  clear  that  in  regions  where  the  solution  is  very  smooth 
(such  as  in  the  shell,  a  reasonable  distance  away  from  the  stiffener),  inserting  a  new 
knot,  and  thus  more  functions  that  maintain  high  continuity,  was  the  more  benefi¬ 
cial  refinement.  In  the  vicinity  of  a  singularity  (such  as  near  the  reentrant  comer 
where  the  shell  meets  the  stiffener  and  the  stress  is  singular),  it  is  more  beneficial 
to  repeat  an  existing  knot  value,  decreasing  the  continuity  of  the  basis  and  simulta¬ 
neously  decreasing  the  support  of  the  basis  functions  in  the  physical  space.  Both  of 
these  effects  help  localize  the  singularity  and  prevent  it  from  polluting  the  results 
elsewhere  in  the  domain 11 . 

The  meshes  for  the  multiple-patch  treatment  of  the  stiffened  shell  are  shown  in 
Figures  44  and  45.  The  locally  refined,  ^-method  meshes  are  seen  in  Figure  45  a.  In 
Figure  45b,  we  see  the  case  where  fewer  elements  are  used.  A  &-type  refinement 
is  used  everywhere  except  at  the  knot  lines  marked  in  red.  The  multiplicities  of 
these  knots  were  increased  with  the  polynomial  order  such  that  the  basis  remained 
C°  across  them.  The  results  for  this  mesh  are  labeled  “Local  &*-ref”  to  indicate 
that  the  ^-refinement  paradigm  was  altered  near  the  singularity.  The  displacements 
are  plotted  versus  the  number  of  degrees  of  freedom  in  Figures  46-49.  The  calcu¬ 
lated  von  Mises  stresses  are  plotted  versus  the  number  of  degrees  of  freedom  in 
Figures  50-53.  The  trunk  space  p-method  results  from  Rank  et  al.  [13]  are  plotted 
for  comparison.  For  displacements,  the  single  patch  results  from  [7]  are  plotted  as 
well. 


4  Conclusions 


We  described  the  possibility  of  h-,  p-,  and  ^-refinement  strategies  and  explored 
their  behavior  on  four  numerical  examples.  The  first  two  examples  concerned  the 
free  vibration  of  an  elastic  rod  and  a  thin  beam.  In  these  cases  the  superiority  of  the 
A  -method  over  the  classical  p-method  is  quite  dramatic.  On  a  per  degree-of-freedom 
basis,  the  results  are  significantly  more  accurate  across  the  entire  spectrum.  In  addi¬ 
tion,  spurious  optical  branches  for  the  p-method  are  found  to  diverge  with  p,  where 
as  optical  branches  are  eliminated  for  the  A:-method.  In  the  latter  case  the  entire 

11  This  is  reminiscent  of  the  heuristic  notion  that  an  hp- method  should  use  large  elements 
with  higher-order  in  smooth  regions  and  small  elements  of  lower-order  near  singularities. 
Coupling  this  with  control  over  the  continuity  across  elements  opens  the  door  to  the  possi¬ 
bility  of  an  hpk- method. 
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Degrees-of-freedom 

Fig.  46.  Hemispherical  shell  with  stiffener.  The  displacement  at  point  A  is  plotted  versus 
the  total  number  of  degrees-of-freedom. 
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Fig.  47.  Hemispherical  shell  with  stiffener.  The  displacement  at  point  B  is  plotted  versus 
the  total  number  of  degrees-of-freedom. 
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Fig.  48.  Hemispherical  shell  with  stiffener.  The  displacement  at  point  C  is  plotted  versus 
the  total  number  of  degrees-of-freedom. 
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Fig.  49.  Hemispherical  shell  with  stiffener.  The  displacement  at  point  D  is  plotted  versus 
the  total  number  of  degrees-of-freedom. 
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Degrees-of-freedom 

Fig.  50.  Hemispherical  shell  with  stiffener.  The  von  Mises  stress  at  point  A  is  plotted  versus 
the  total  number  of  degrees-of-freedom. 
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Fig.  51.  Hemispherical  shell  with  stiffener.  The  von  Mises  stress  at  point  B  is  plotted  versus 
the  total  number  of  degrees-of-freedom. 
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Fig.  52.  Hemispherical  shell  with  stiffener.  The  von  Mises  stress  at  point  C  is  plotted  versus 
the  total  number  of  degrees-of-freedom. 
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Fig.  53.  Hemispherical  shell  with  stiffener.  The  von  Mises  stress  at  point  D  is  plotted  versus 
the  total  number  of  degrees-of-freedom. 
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spectrum  converges  as  order  is  increased.  We  feel  in  these  cases  that  the  increased 
smoothness  of  the  /.-method  basis  functions  better  exploits  the  smoothness  of  the 
exact  analytical  eigensolutions  than  do  the  C'°-continuous  basis  functions  of  the  p- 
method.  For  further  studies  of  the  behavior  of  isogeometric  approaches  to  structural 
vibrations,  see  Cottrell  et  al.  [3]. 

We  then  considered  two  elliptic  boundary  value  problems  for  shell  structures  ex¬ 
hibiting  singular  behavior.  The  first  shell  was  a  hyperboloid  subjected  to  circum¬ 
ferentially  varying  pressure  (Lee  and  Bathe  [10]).  There  is  a  weak  boundary  layer 
singularity.  The  second  shell  was  a  spherical  cap,  with  a  hole  at  the  apex,  built  into 
a  solid  ring  stiffener  (Rank  et  al.  [13]).  The  intersection  of  the  shell  and  the  stiffener 
creates  a  strong  line  singularity.  These  problems  provided  the  opportunity  to  study 
the  effects  of  smoothness  of  basis  functions  in  the  vicinity  of  singularities.  In  both 
cases,  we  used  a  trivariate  NURBS  solid  description  instead  of  a  shell  theory.  For 
the  hyperboloidal  shell,  we  reduced  smoothness  locally  in  the  vicinity  of  the  bound¬ 
ary  layer.  In  the  case  of  coarse  meshes,  this  improved  accuracy,  where  as  for  finer 
meshes  the  pure  A;-method  was  more  accurate.  In  the  case  of  the  spherical  cap,  we 
employed  a  multipatch  approach  with  local  refinement  and  again  compared  smooth 
discretizations  within  the  patches  with  ones  in  which  continuity  was  reduced  to  C° 
in  the  vicinity  of  the  singularity.  We  found  this  latter  approach  led  to  more  rapid 
convergence.  The  reason  for  this  seems  to  be  that  basis  functions  having  support  in 
the  vicinity  of  the  singularity  tend  to  propagate  information  away  from  the  singu¬ 
larity.  The  support  of  smooth  A;-method  basis  functions  is  greater  than  the  support 
of  the  same  order  p-method  functions  when  there  are  approximately  the  same  num¬ 
bers  of  degrees-of-freedom.  As  a  result,  the  errors  created  by  the  singularities  tend 
to  propagate  further  for  the  smoother  basis  functions  of  the  A;-method.  By  judi¬ 
ciously  locating  a  few  surfaces  of  reduced  continuity,  the  “pollution”  created  by 
the  singularities  seemed  to  be  more  locally  confined. 

In  conclusion,  or  studies  revealed  some  cases  where  increased  smoothness  dramat¬ 
ically  improved  accuracy  and  other  cases  where  some  local  reduction  in  smooth¬ 
ness  also  enhanced  accuracy.  Historically,  the  dominance  of  C°  finite  elements  has 
precluded  the  opportunity  for  assessing  the  potential  of  increased  smoothness.  It 
is  clear  from  the  studies  described  herein,  and  those  in  our  other  recent  papers 
[1,  3,  7],  that  the  potential  is  very  significant.  However,  it  does  not  seem  to  be  a 
black  and  white  issue,  but  rather  to  depend  strongly  on  the  nature  of  the  exact  solu¬ 
tion.  Further  studies  need  to  be  performed  to  assess  the  tradeoffs  in  a  wider  variety 
of  problems. 
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Richardson  Extrapolation 
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Richardson  extrapolation  is  a  technique  in  which  two  approximations  of  a  known 
order  of  accuracy  relative  to  some  parameter  are  combined  to  obtain  an  approxima¬ 
tion  with  a  higher  order  of  accuracy.  In  the  case  of  finite  elements,  the  parameter  is 
the  mesh  size,  h,  and  the  order  of  accuracy  is  dictated  by  the  polynomial  order  of 
the  basis  functions 12  .  In  general,  we  generate  an  approximation  A(h)  to  a  desired 
quantity  A,  where  the  order  of  the  error  is  known.  That  is, 

A  =  A(h)  +  Bhk  +  Chk+1  +  Dhk+2  +  ...  (A.l) 

where  k  is  known  but  B ,  C.  D.  etc.  are  unknown  constants.  More  concisely,  we 
write 

A  =  A(h)  +  Bhk  +  0(hk+1).  (A.2) 

We  can  get  a  second  such  equation  by  generating  a  second  approximation  with  a 
different  mesh  size,  for  example 

A  =  A(h/2)  +  B(h/2)k +  0(hk+1).  (A.3) 

We  can  remove  the  lowest  order  term  in  the  error  by  taking  2k  times  A.3  and  sub¬ 
tracting  A.2,  which  yields 

(: 2k  -  1  )A  =  2kA(h/2)  -  A(h)  +  0(hk+1).  (A.4) 


Simplifying,  we  arrive  at 

A  =  —  +  0(hk+l). 


(A.5) 


We  now  have  a  new  approximation 

a  .  (A.6) 

whose  order  of  accuracy  is  higher  than  either  of  the  two  approximations  that  gen¬ 
erated  it. 
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